Aakriti Poudel
  • Home
  • About
  • Resources
  • Posts

On this page

  • Load required libraries
  • Read data
  • Merge raster data and create houston bounding box
    • Visualize before and after night light intensities
  • Create blackout mask
  • Estimate of the number of homes in Houston that lost power
    • Map of the homes in Houston that lost power
  • Census tracts in Houston that lost power
    • Map of the census tracts in Houston that lost power
  • Plot comparing the distributions of median household income
  • Brief reflection

Identifying the impacts of extreme weather

MEDS
R
GIS
geospatial
remote sensing
data science
A case study of Houston, Texas
Author

Aakriti Poudel

Published

December 3, 2025

Load required libraries

Code
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.3.0
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.1.0     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code
library(sf)
Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
Code
library(tmap)
library(terra)
terra 1.8.60

Attaching package: 'terra'

The following object is masked from 'package:tidyr':

    extract
Code
library(stars)
Loading required package: abind

Read data

Code
# Read night lights data
tile05_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v05.001.2021039064328.tif'))
tile06_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v06.001.2021039064329.tif'))
tile05_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather','data', 'VNP46A1', 'VNP46A1.A2021047.h08v05.001.2021048091106.tif'))
tile06_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021047.h08v06.001.2021048091105.tif'))


# Read roads data
highways <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'gis_osm_roads_free_1.gpkg'),
                 query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass = 'motorway'") %>%
  st_transform(crs = 'EPSG:3083')


# Read houses data
houses <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data',
                          'gis_osm_buildings_a_free_1.gpkg'),
  query = "SELECT *
    FROM gis_osm_buildings_a_free_1
    WHERE (type IS NULL AND name IS NULL)
      OR type IN ('residential', 'apartments', 'house', 'static_caravan', 'detached')") %>% 
    st_transform('EPSG:3083')


# Explore the contents of the geodatabase
socioeconomic <- st_layers(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'))

# Extract the census tract information
census_tract <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'ACS_2019_5YR_TRACT_48_TEXAS') %>% 
  st_transform('EPSG:3083')

# Extract the income information
income <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'X19_INCOME')

Merge raster data and create houston bounding box

Code
# Merge two raster data of February 7
nightlight_07 <- st_mosaic(tile05_07, tile06_07)

# Merge two raster data of February 16
nightlight_16 <- st_mosaic(tile05_16, tile06_16)

# Create Houston bounding box
houston_bbox <- st_bbox(c(xmin = -96.5, xmax = -94.5, ymin = 29, ymax = 30.5),
                        crs = 'EPSG:4326')

# Crop all raster data to Houston area
nightlight_07_crop <- st_crop(nightlight_07, houston_bbox)
nightlight_07_crop[(nightlight_07_crop > 1000) | (nightlight_07_crop <= 0)] <- NA

nightlight_16_crop <- st_crop(nightlight_16, houston_bbox)
nightlight_16_crop[(nightlight_16_crop > 1000) | (nightlight_16_crop <= 0)] <- NA

Visualize before and after night light intensities

Code
# Map for February 7 (before storm)
beforestorm <- tm_shape(nightlight_07_crop) +
  tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
            col.legend = tm_legend(title = 'Night light\nintensity',
                                   position = tm_pos_out('right'))) +
  tm_title('Houston Night Lights\nFebruary 07, 2021 (Before storm)',
           size = 1.3) +
  tm_compass(position = c('left', 'top'),
             type = 'arrow',
             size = 3,
             text.size = 0.6,
             color.dark = 'grey50',
             text.color = 'white') +  
  tm_scalebar(position = c('right', 'bottom'),
              text.size = 0.8,
              color.dark = 'grey50',
              text.color = 'white') +  
  tm_layout(bg.color = 'white',           
            outer.bg.color = 'white',      
            frame = FALSE,
            legend.show = FALSE)


# Map for February 16 (after storm)
afterstorm <- tm_shape(nightlight_16_crop) +
  tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
            col.legend = tm_legend(title = 'Night light\nintensity',
                                   position = tm_pos_out('right'))) +
  tm_title('Houston Night Lights\nFebruary 16, 2021 (After storm)',
           size = 1.3) +
  tm_compass(position = c('left', 'top'),
             type = 'arrow',
             size = 3,
             text.size = 0.6,
             color.dark = 'grey50',
             text.color = 'white') +  
  tm_scalebar(position = c('right', 'bottom'),
              text.size = 0.8,
              color.dark = 'grey50',
              text.color = 'white') +  
  tm_layout(frame = FALSE,
            legend.show = FALSE)

legend_map <- tm_shape(nightlight_16_crop) +
  tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
            col.legend = tm_legend(title = 'Night light\nintensity',
                                   position = tm_pos_out('right'))) +
  tm_layout(legend.only = TRUE, legend.outside = TRUE, legend.position = c("right", "center"))

tmap_arrange(beforestorm, afterstorm, legend_map, ncol = 3)

Create blackout mask

Code
# Calculate the difference in light intensity
nightlight_diff <- nightlight_07 - nightlight_16

# Crop and reclassify the raster data 
mask <- st_crop(nightlight_diff, houston_bbox)
mask[mask < 200] <- NA

# Vectorize the blackout mask
blackout <- st_as_sf(mask,
                     as_points = FALSE,
                     merge = TRUE) %>%
  st_make_valid() %>%
  st_transform(crs = 'EPSG:3083')
Code
# Combine all highway geometries into one
highways_union <- st_union(highways)

# Create 200m buffer around all highways
highways_buffer <- st_buffer(highways_union, dist = 200)

# Find areas that experienced blackouts further than 200m from highways
blackout_far <- st_difference(blackout, highways_buffer)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Code
# Visualize the blackout in Houston
tm_basemap('OpenStreetMap') +
  tm_shape(blackout_far) +
    tm_polygons(fill = 'ivory', fill_alpha = 0.8, col = 'black') +
    tm_title('Blackout areas in Houston', fontface = 'bold', size = 1) +
    tm_layout(legend.show = TRUE) +
    tm_scalebar(position = c('left', 'bottom')) +
    tm_compass(position = c('right', 'top'))

Estimate of the number of homes in Houston that lost power

Approximately 157,970 homes in Houston lost power in February 2021 due to the severe winter storms and resulting major power crisis.

Note: I used st_intersects. st_intersects returns a boolean value indicating whether two geometries share any features or if they have any point in common. TRUE (or 1) if the geometries intersect, FALSE (or 0) if they are completely separate.

Code
# Keep buildings that intersect with blackout areas
houses_blackout <- st_intersects(houses, blackout_far)

# Check if each building is in a blackout area
in_blackout <- lengths(houses_blackout) > 0

# Filter buildings in blackout areas
blackout_houses <- houses[in_blackout, ]

# Number of impacted buildings
n_blackout_houses <- nrow(blackout_houses)

# Convert all buildings to centroids
houses_points <- houses %>% 
  st_centroid()
Warning: st_centroid assumes attributes are constant over geometries
Code
# Houses impacted by blackouts
houses_impacted <- blackout_houses %>% 
  st_centroid()
Warning: st_centroid assumes attributes are constant over geometries
Code
# Houses not impacted by blackouts
houses_not_impacted <- houses_points %>% 
  filter(!osm_id %in% blackout_houses$osm_id)

Map of the homes in Houston that lost power

Code
# Create the map
tm_basemap('OpenStreetMap') +
  tm_shape(houses_not_impacted) +
    tm_dots(fill = 'midnightblue', size = 0.02, fill_alpha = 0.8) +
  tm_shape(houses_impacted) +
    tm_dots(fill = 'firebrick', size = 0.02, fill_alpha = 0.8) + 
  tm_title('Houses in Houston that lost power',
           position = tm_pos_out('center', 'top'), 
           fontface = 'bold', size = 1.2) +
  tm_layout(legend.outside = TRUE,
            legend.outside.position = 'right',
            legend.title.size = 1,
            legend.text.size = 0.8) +
  tm_add_legend(labels = c('Impacted houses', 'Not impacted houses'),
                fill = c('firebrick', 'midnightblue'),
                col = c('firebrick', 'midnightblue'),
                alpha = c(0.2, 0.7, 0.5),
                title = 'Legend') +
  tm_scalebar(position = c('left', 'bottom'),
              text.size = 0.6) +
  tm_compass(position = c('right', 'bottom'),
             type = '4star',
             size = 2)

Census tracts in Houston that lost power

Code
# Rename the column in income dataset to match with the census tract
income_renamed <- income %>% 
  rename("GEOID_Data" = "GEOID")

# Left join the two datasets by the geometry
census_income <- census_tract %>% left_join(income_renamed, by = 'GEOID_Data')
Code
# Check CRS of the dataset
if(st_crs(census_income) == st_crs(houston_bbox)){
print("Coordinate reference systems match!")
} else {
warning("Update coordinate reference systems to match!")
}
Warning: Update coordinate reference systems to match!
Code
# Transform the CRS and crop to houston bounding box
census_income <- st_transform(census_income, crs = 'epsg:4326') 

# Crop census income to Houston
census_income_crop <- census_income %>% 
  st_crop(houston_bbox)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Code
# Transform and clean blackout houses
blackout_houses_transformed <- blackout_houses %>% 
  st_transform(crs = 'epsg:4326')

# Identify the blackout areas
census_blackout <- st_intersects(census_income_crop, blackout_houses_transformed)

# Create the column
census_blackout_col <- census_income_crop %>% 
  mutate(blackout = lengths(census_blackout) > 0)

# Select only the values that are true for plotting
census_blackout_col_filtered <- census_blackout_col[census_blackout_col$blackout, ]

Map of the census tracts in Houston that lost power

Code
# Map the census tracts that lost power
tm_basemap('CartoDB.VoyagerNoLabels') +
  tm_shape(census_blackout_col_filtered) +
  tm_polygons(fill = 'blackout',
              fill.scale = tm_scale_categorical(
                values = c('TRUE' = 'firebrick'),
                labels = c('TRUE' = 'Blackout')),
              fill.legend = tm_legend(title = 'Blackout status'),
              col = 'grey50',
              lwd = 0.2) +
  tm_title('Census tracts in Houston that lost power',
           position = tm_pos_out('center', 'top'),
           fontface = 'bold',
           size = 1) +
  tm_layout(legend.outside = TRUE,
            legend.outside.position = 'right',
            legend.bg.color = 'white',
            legend.bg.alpha = 0.8) +
  tm_compass(position = c('right', 'bottom'),
             type = 'arrow',
             size = 1) +
  tm_scalebar(position = c('left', 'top'),
              text.size = 0.8)

Plot comparing the distributions of median household income

Code
ggplot(census_blackout_col %>% filter(!is.na(B19013e1)), # Remove NA values
                                      aes( x = blackout, y = B19013e1, fill = blackout)) +
  geom_boxplot(alpha = 0.8) +
  scale_x_discrete(labels = c('TRUE' = 'Blackout',
                              'FALSE' = 'No blackout')) +
  scale_fill_manual(values = c('TRUE' = 'firebrick',
                               'FALSE' = 'midnightblue'),
                    labels = c('Blackout', 'No blackout')) +
  labs(title = 'Median household income for census tracts',
       subtitle = 'Blackout status in Houston census tracts',
       x = 'Blackout status',
       y = 'Median household income in $',
       fill = 'Status') +
  scale_y_continuous(labels = scales::dollar_format()) +
  theme_minimal() +
  theme(plot.title = element_text(face = 'bold', size = 16),
        plot.subtitle = element_text(size = 14),
        legend.position = 'bottom')
Figure 1: This plot compares the distributions of median household income between census tracts that experienced blackouts and those that did not during the severe winter storms of February 2021 in Houston, Texas.

Brief reflection

A brief reflection (approx. 100 words) summarizing your results and discussing any limitations to this study

Answer This analysis used satellite data to compare median household incomes between areas affected and unaffected by the blackout. Interestingly, higher-income census tracts were more affected, as their median income was slightly higher than that of non-blackout areas. This result may reflect several limitations such as wealthier neighborhoods often have more outdoor lighting, census tract–level data may overlook neighborhood differences, and blackout areas could be misclassified due to highway lighting or cloud cover. The analysis may also have assumed all structures were residential, while vegetation near power lines and the 200-meter buffer used could further limit accuracy of true lighting.

Citation

BibTeX citation:
@online{poudel2025,
  author = {Poudel, Aakriti},
  title = {Identifying the Impacts of Extreme Weather},
  date = {2025-12-03},
  url = {https://aakriti-poudel-chhetri.github.io/posts/2025-12-impacts-of-extreme-weather/},
  langid = {en}
}
For attribution, please cite this work as:
Poudel, Aakriti. 2025. “Identifying the Impacts of Extreme Weather.” December 3, 2025. https://aakriti-poudel-chhetri.github.io/posts/2025-12-impacts-of-extreme-weather/.

© 2025, Aakriti Poudel

 

This website is built with Quarto, R and GitHub.